Chapter 6 Diversity analysis

6.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db2) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

6.1.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.3 Post-Transplant_1 samples

bxp1<-alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
        scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")


# Add significance bars with p-values
bxp1 + 
  geom_signif(comparisons = list(c("Control", "Hot_control"), c("Control", "Treatment"), c("Hot_control", "Treatment")),
              map_signif_level = TRUE)

6.1.3.1 Richness

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point =="5_Post-FMT1") 


Modelq0GLMMNB <- glmer.nb(richness ~ type+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(15311.33)  ( log )
Formula: richness ~ type + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   248.5    254.8   -119.3    238.5       21 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.35323 -0.19624  0.00957  0.16293  0.31217 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.2757   0.5251  
Number of obs: 26, groups:  individual, 26

Fixed effects:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)       3.5033     0.1854  18.893   <2e-16 ***
typeHot_control   0.5046     0.2590   1.949   0.0513 .  
typeTreatment     0.2199     0.2692   0.817   0.4140    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) typHt_
typHt_cntrl -0.716       
typeTretmnt -0.686  0.491

6.1.3.2 Neutral

Model_neutral <- lme(fixed = neutral ~ type, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  188.7488 194.4263 -89.37439

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    9.586725 3.595032

Fixed effects:  neutral ~ type 
                    Value Std.Error DF  t-value p-value
(Intercept)     17.359333  3.412877 23 5.086423  0.0000
typeHot_control  8.380488  4.826537 23 1.736336  0.0959
typeTreatment    7.504392  4.975080 23 1.508396  0.1451
 Correlation: 
                (Intr) typHt_
typeHot_control -0.707       
typeTreatment   -0.686  0.485

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-0.70984520 -0.23549305 -0.03481478  0.24385684  0.54281137 

Number of Observations: 26
Number of Groups: 26 

6.1.3.3 Phylogenetic

Model_phylo <- lme(fixed = phylogenetic ~ Population, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  72.19064 76.90286 -32.09532

Random effects:
 Formula: ~1 | individual
        (Intercept)  Residual
StdDev:   0.7770718 0.2914019

Fixed effects:  phylogenetic ~ Population 
                      Value Std.Error DF  t-value p-value
(Intercept)        4.307539 0.2012835 24 21.40036  0.0000
PopulationHot_dry -0.036754 0.3421161 24 -0.10743  0.9153
 Correlation: 
                  (Intr)
PopulationHot_dry -0.588

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-0.71953431 -0.17146063  0.01286462  0.19655255  0.64008560 

Number of Observations: 26
Number of Groups: 26 

6.1.3.4 Functional

Model_func <- lme(fixed = functional ~ type, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_func)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
        AIC       BIC   logLik
  -27.23436 -21.55689 18.61718

Random effects:
 Formula: ~1 | individual
        (Intercept)   Residual
StdDev:  0.08760598 0.03285224

Fixed effects:  functional ~ type 
                    Value  Std.Error DF  t-value p-value
(Intercept)     1.4854053 0.03118774 23 47.62786  0.0000
typeHot_control 0.0073281 0.04410613 23  0.16615  0.8695
typeTreatment   0.0032476 0.04546356 23  0.07143  0.9437
 Correlation: 
                (Intr) typHt_
typeHot_control -0.707       
typeTreatment   -0.686  0.485

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-0.73243259 -0.23582828 -0.03001854  0.23261622  0.61672731 

Number of Observations: 26
Number of Groups: 26 

6.1.4 Post-Transplant_2 samples

bxp1<-alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
          scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

# Add significance bars with p-values
bxp1 + 
  geom_signif(comparisons = list(c("Control", "Hot_control"), c("Control", "Treatment"), c("Hot_control", "Treatment")),
              map_signif_level = TRUE)

6.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

6.3 Permanovas

6.3.0.1 Load required data

6.3.1 1. Are the wild populations similar?

6.3.1.1 Wild: cold vs warm adapted lizards

6.3.1.1.1 Number of samples used

[1] 27

6.3.1.2 Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000012 0.000012 0.0012    999  0.978
Residuals 25 0.257281 0.010291                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.98
Hot_dry   0.97302        

adonis2(richness ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.542719	0.2095041	6.625717	0.001
Residual	25	5.820951	0.7904959	NA	NA
Total	26	7.363669	1.0000000	NA	NA

6.3.1.3 Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000048 0.0000476 0.0044    999  0.951
Residuals 25 0.270114 0.0108046                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.954
Hot_dry   0.94763        

adonis2(neutral ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.918266	0.2608511	8.822682	0.001
Residual	25	5.435610	0.7391489	NA	NA
Total	26	7.353876	1.0000000	NA	NA

6.3.1.4 Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03585 0.035847 2.4912    999  0.133
Residuals 25 0.35973 0.014389                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.126
Hot_dry   0.12705        

adonis2(phylo ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.3218613	0.2162815	6.899207	0.001
Residual	25	1.1662981	0.7837185	NA	NA
Total	26	1.4881594	1.0000000	NA	NA

6.3.1.5 Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.018367 0.018367 1.5597    999  0.209
Residuals 25 0.294402 0.011776                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.224
Hot_dry   0.22328        

adonis2(func ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(func))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.08585779	0.172879	5.225323	0.063
Residual	25	0.41077745	0.827121	NA	NA
Total	26	0.49663525	1.000000	NA	NA

beta_q0n_nmds_wild <- richness %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta, by = join_by(sample == Tube_code))

beta_q1n_nmds_wild <- neutral %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta, by = join_by(sample == Tube_code))

beta_q1p_nmds_wild <- phylo %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta, by = join_by(sample == Tube_code))

beta_q1f_nmds_wild <- func %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta, by = join_by(sample == Tube_code))

6.3.2 2. Effect of acclimation

6.3.2.0.1 Number of samples used

[1] 26

6.3.2.1 Richness

richness_accli <- as.matrix(beta_q0n$S)
richness_accli <- as.dist(richness_accli[rownames(richness_accli) %in% samples_to_keep_accli,
               colnames(richness_accli) %in% samples_to_keep_accli])
betadisper(richness_accli, subset_meta_accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.093187 0.093187 11.812    999  0.004 **
Residuals 24 0.189340 0.007889                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.002
Hot_dry  0.0021532        

adonis2(richness_accli ~ Population,
        data = subset_meta_accli %>% arrange(match(Tube_code,labels(richness_accli))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.639630	0.1929088	5.736415	0.001
Residual	24	6.859879	0.8070912	NA	NA
Total	25	8.499509	1.0000000	NA	NA

6.3.2.2 Neutral

neutral_accli <- as.matrix(beta_q1n$S)
neutral_accli <- as.dist(neutral_accli[rownames(neutral_accli) %in% samples_to_keep_accli,
               colnames(neutral_accli) %in% samples_to_keep_accli])
betadisper(neutral_accli, subset_meta_accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.05603 0.056026 4.1918    999  0.041 *
Residuals 24 0.32077 0.013365                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.043
Hot_dry  0.051717        

adonis2(neutral_accli ~ Population,
        data = subset_meta_accli %>% arrange(match(Tube_code,labels(neutral_accli))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.961192	0.2493171	7.970889	0.001
Residual	24	5.905063	0.7506829	NA	NA
Total	25	7.866255	1.0000000	NA	NA

6.3.2.3 Phylogenetic

phylo_accli <- as.matrix(beta_q1p$S)
phylo_accli <- as.dist(phylo_accli[rownames(phylo_accli) %in% samples_to_keep_accli,
               colnames(phylo_accli) %in% samples_to_keep_accli])
betadisper(phylo_accli, subset_meta_accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03637 0.036365 2.3087    999  0.142
Residuals 24 0.37804 0.015752                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.15
Hot_dry   0.14172        

adonis2(phylo_accli ~ Population,
        data = subset_meta_accli %>% arrange(match(Tube_code,labels(phylo_accli))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.2113395	0.1379515	3.84066	0.009
Residual	24	1.3206449	0.8620485	NA	NA
Total	25	1.5319844	1.0000000	NA	NA

6.3.2.4 Functional

func_accli <- as.matrix(beta_q1f$S)
func_accli <- as.dist(func_accli[rownames(func_accli) %in% samples_to_keep_accli,
               colnames(func_accli) %in% samples_to_keep_accli])
betadisper(func_accli, subset_meta_accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00471 0.0047108 0.2042    999  0.662
Residuals 24 0.55370 0.0230708                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.664
Hot_dry   0.65542        

adonis2(func_accli ~ Population,
        data = subset_meta_accli %>% arrange(match(Tube_code,labels(func_accli))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.00470883	0.006952167	0.1680201	0.614
Residual	24	0.67260942	0.993047833	NA	NA
Total	25	0.67731825	1.000000000	NA	NA

beta_q0n_nmds_accli <- richness_accli %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_accli, by = join_by(sample == Tube_code))

beta_q1n_nmds_accli <- neutral_accli %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_accli, by = join_by(sample == Tube_code))

beta_q1p_nmds_accli <- phylo_accli %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_accli, by = join_by(sample == Tube_code))

beta_q1f_nmds_accli <- func_accli %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_accli, by = join_by(sample == Tube_code))

6.3.3 3. Comparison between Wild and Acclimation

6.3.3.0.1 Number of samples used

[1] 53

6.3.3.0.2 Richness

richness_accli1 <- as.matrix(beta_q0n$S)
richness_accli1 <- as.dist(richness_accli1[rownames(richness_accli1) %in% samples_to_keep_accli1,
               colnames(richness_accli1) %in% samples_to_keep_accli1])
betadisper(richness_accli1, subset_meta_accli1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.03286 0.032865 2.9698    999  0.094 .
Residuals 51 0.56438 0.011066                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.083
Hot_dry   0.09089        

adonis2(richness_accli1 ~ Population*time_point,
        data = subset_meta_accli1 %>% arrange(match(Tube_code,labels(richness_accli1))),
        permutations = 999,
        strata = subset_meta_accli1 %>% arrange(match(Tube_code,labels(richness_accli1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	3.762933	0.2288365	4.846784	0.001
Residual	49	12.680830	0.7711635	NA	NA
Total	52	16.443763	1.0000000	NA	NA

#Arrange of metadata dataframe
subset_meta_accli1_arrange <- column_to_rownames(subset_meta_accli1, "Tube_code")
subset_meta_accli1_arrange<-subset_meta_accli1_arrange[labels(richness_accli1),]

pairwise<-pairwise.adonis(richness_accli1,subset_meta_accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation	1	1.6396299	5.736415	0.19290877	0.001	0.006	*
Cold_wet.1_Acclimation vs Cold_wet.0_Wild	1	0.5273732	1.906862	0.05462715	0.010	0.060
Cold_wet.1_Acclimation vs Hot_dry.0_Wild	1	1.5558412	5.190707	0.17782052	0.001	0.006	*
Hot_dry.1_Acclimation vs Cold_wet.0_Wild	1	1.8259388	8.319131	0.24968031	0.001	0.006	*
Hot_dry.1_Acclimation vs Hot_dry.0_Wild	1	0.3856333	1.736034	0.09788177	0.009	0.054
Cold_wet.0_Wild vs Hot_dry.0_Wild	1	1.5427188	6.625717	0.20950408	0.001	0.006	*

6.3.3.0.3 Neutral

neutral_accli1 <- as.matrix(beta_q1n$S)
neutral_accli1 <- as.dist(neutral_accli1[rownames(neutral_accli1) %in% samples_to_keep_accli1,
               colnames(neutral_accli1) %in% samples_to_keep_accli1])
betadisper(neutral_accli1, subset_meta_accli1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01431 0.014310 1.1177    999  0.313
Residuals 51 0.65296 0.012803                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.297
Hot_dry   0.29539        

adonis2(neutral_accli1 ~ Population*time_point,
        data = subset_meta_accli1 %>% arrange(match(Tube_code,labels(neutral_accli1))),
        permutations = 999,
        strata = subset_meta_accli1 %>% arrange(match(Tube_code,labels(neutral_accli1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	4.735485	0.2945657	6.820252	0.001
Residual	49	11.340673	0.7054343	NA	NA
Total	52	16.076158	1.0000000	NA	NA

pairwise<-pairwise.adonis(neutral_accli1,subset_meta_accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation	1	1.9611917	7.970889	0.24931708	0.001	0.006	*
Cold_wet.1_Acclimation vs Cold_wet.0_Wild	1	0.7466371	3.064914	0.08498327	0.001	0.006	*
Cold_wet.1_Acclimation vs Hot_dry.0_Wild	1	2.1421077	8.176866	0.25412251	0.001	0.006	*
Hot_dry.1_Acclimation vs Cold_wet.0_Wild	1	2.0371666	10.078295	0.28730857	0.001	0.006	*
Hot_dry.1_Acclimation vs Hot_dry.0_Wild	1	0.6314393	3.060027	0.16054681	0.001	0.006	*
Cold_wet.0_Wild vs Hot_dry.0_Wild	1	1.9182663	8.822682	0.26085105	0.001	0.006	*

6.3.3.0.4 Phylogenetic

phylo_accli1 <- as.matrix(beta_q1p$S)
phylo_accli1 <- as.dist(phylo_accli1[rownames(phylo_accli1) %in% samples_to_keep_accli1,
               colnames(phylo_accli1) %in% samples_to_keep_accli1])
betadisper(phylo_accli1, subset_meta_accli1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     1 0.00001 0.0000111 6e-04    999  0.984
Residuals 51 0.89017 0.0174543                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.985
Hot_dry   0.97998        

adonis2(phylo_accli1 ~ Population*time_point,
        data = subset_meta_accli1 %>% arrange(match(Tube_code,labels(phylo_accli1))),
        permutations = 999,
        strata = subset_meta_accli1 %>% arrange(match(Tube_code,labels(phylo_accli1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	0.7622566	0.2345983	5.006223	0.001
Residual	49	2.4869430	0.7654017	NA	NA
Total	52	3.2491997	1.0000000	NA	NA

pairwise<-pairwise.adonis(phylo_accli1,subset_meta_accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation	1	0.2113395	3.840660	0.1379515	0.010	0.060
Cold_wet.1_Acclimation vs Cold_wet.0_Wild	1	0.2076830	4.216762	0.1133028	0.003	0.018	.
Cold_wet.1_Acclimation vs Hot_dry.0_Wild	1	0.3635093	4.799786	0.1666605	0.001	0.006	*
Hot_dry.1_Acclimation vs Cold_wet.0_Wild	1	0.2121479	7.924059	0.2406769	0.002	0.012	.
Hot_dry.1_Acclimation vs Hot_dry.0_Wild	1	0.2092433	3.885515	0.1953942	0.001	0.006	*
Cold_wet.0_Wild vs Hot_dry.0_Wild	1	0.3218613	6.899207	0.2162815	0.001	0.006	*

6.3.3.0.5 Functional

func_accli1 <- as.matrix(beta_q1f$S)
func_accli1 <- as.dist(func_accli1[rownames(func_accli1) %in% samples_to_keep_accli1,
               colnames(func_accli1) %in% samples_to_keep_accli1])
betadisper(func_accli1, subset_meta_accli1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01354 0.013538 0.7762    999   0.41
Residuals 51 0.88953 0.017442                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.408
Hot_dry   0.38245        

adonis2(func_accli1 ~ Population*time_point,
        data = subset_meta_accli1 %>% arrange(match(Tube_code,labels(func_accli1))),
        permutations = 999,
        strata = subset_meta_accli1 %>% arrange(match(Tube_code,labels(func_accli1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	0.1059556	0.08908753	1.597405	0.206
Residual	49	1.0833869	0.91091247	NA	NA
Total	52	1.1893425	1.00000000	NA	NA

pairwise<-pairwise.adonis(func_accli1,subset_meta_accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation	1	0.004708830	0.16802012	0.006952167	0.603	1.000
Cold_wet.1_Acclimation vs Cold_wet.0_Wild	1	0.001001660	0.03572249	0.001081329	0.691	1.000
Cold_wet.1_Acclimation vs Hot_dry.0_Wild	1	0.088628476	3.78330596	0.136171914	0.092	0.552
Hot_dry.1_Acclimation vs Cold_wet.0_Wild	1	0.002762199	0.13250295	0.005272175	0.605	1.000
Hot_dry.1_Acclimation vs Hot_dry.0_Wild	1	0.042282101	4.27996390	0.211043961	0.096	0.576
Cold_wet.0_Wild vs Hot_dry.0_Wild	1	0.085857794	5.22532292	0.172878978	0.055	0.330

beta_richness_nmds_accli1 <- richness_accli1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_accli1, by = c("sample" = "Tube_code"))

beta_neutral_nmds_accli1 <- neutral_accli1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_accli1, by = c("sample" = "Tube_code"))

beta_phylo_nmds_accli1 <- phylo_accli1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_accli1, by = join_by(sample == Tube_code))

beta_func_nmds_accli1 <- func_accli1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_accli1, by = join_by(sample == Tube_code))

6.3.4 4. Effect of FMT on microbiota community

6.3.4.1 Comparison between Acclimation vs Post-FMT1

6.3.4.1.1 Number of samples used

[1] 52

6.3.4.1.2 Richness

richness_post3 <- as.matrix(beta_q0n$S)
richness_post3 <- as.dist(richness_post3[rownames(richness_post3) %in% samples_to_keep_post3,
               colnames(richness_post3) %in% samples_to_keep_post3])
betadisper(richness_post3, subset_meta_post3$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.10843 0.108427 25.578    999  0.001 ***
Residuals 50 0.21195 0.004239                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Cold_wet Hot_dry
Cold_wet              0.001
Hot_dry  6.0928e-06        

adonis2(richness_post3 ~ Population*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	3.472191	0.193401	3.836375	0.001
Residual	48	14.481131	0.806599	NA	NA
Total	51	17.953321	1.000000	NA	NA

#Arrange of metadata dataframe
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(richness_post3),]

pairwise <- pairwise.adonis(richness_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.3657243	1.123239	0.06966584	0.258	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.2605860	4.943717	0.24788342	0.001	0.015	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.6804988	2.100611	0.12283837	0.001	0.015	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8481258	2.673290	0.16033371	0.001	0.015	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.1180661	3.809275	0.20252111	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.3606630	5.087152	0.24124415	0.001	0.015	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.7216200	2.172734	0.11956009	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9551308	2.926054	0.16322910	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.2263345	4.039487	0.20157637	0.001	0.015	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.4319792	5.384836	0.25180628	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8172413	3.194690	0.17558364	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.5796135	2.441615	0.13239702	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.5615418	1.729004	0.10335366	0.013	0.195
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.8438429	2.793772	0.14865413	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.3734921	1.268929	0.07799710	0.108	1.000

6.3.4.1.3 Neutral

neutral_post3 <- as.matrix(beta_q1n$S)
neutral_post3 <- as.dist(neutral_post3[rownames(neutral_post3) %in% samples_to_keep_post3,
               colnames(neutral_post3) %in% samples_to_keep_post3])
betadisper(neutral_post3, subset_meta_post3$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.09089 0.090889 12.898    999  0.002 **
Residuals 50 0.35233 0.007047                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Cold_wet Hot_dry
Cold_wet              0.001
Hot_dry  0.00074926        

adonis2(neutral_post3 ~ Population*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	4.443476	0.2628643	5.705637	0.001
Residual	48	12.460591	0.7371357	NA	NA
Total	51	16.904067	1.0000000	NA	NA

pairwise <- pairwise.adonis(neutral_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.2759858	0.9928976	0.06208366	0.457	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.4146464	6.5078325	0.30257965	0.001	0.015	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.8153894	3.0970603	0.17113610	0.001	0.015	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	1.1809241	4.4856470	0.24265567	0.001	0.015	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.4321524	5.7774260	0.27806264	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.6347704	6.8326887	0.29925029	0.001	0.015	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.9517634	3.3715700	0.17404733	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	1.3127773	4.6298256	0.23585668	0.002	0.030	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.6713369	6.2395460	0.28056085	0.001	0.015	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.5409781	6.8338056	0.29928456	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9133614	4.0964534	0.21451383	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.6954835	3.2951234	0.17077493	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.6051778	2.2508491	0.13047758	0.011	0.165
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	1.0528902	4.1436369	0.20570451	0.002	0.030	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.4150076	1.6372683	0.09840968	0.049	0.735

6.3.4.1.4 Phylogenetic

phylo_post3 <- as.matrix(beta_q1p$S)
phylo_post3 <- as.dist(phylo_post3[rownames(phylo_post3) %in% samples_to_keep_post3,
               colnames(phylo_post3) %in% samples_to_keep_post3])
betadisper(phylo_post3, subset_meta_post3$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.04088 0.040882 2.9254    999  0.097 .
Residuals 50 0.69874 0.013975                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.099
Hot_dry  0.093394        

adonis2(phylo_post3 ~ Population*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	0.6468577	0.2180259	4.461035	0.001
Residual	48	2.3200272	0.7819741	NA	NA
Total	51	2.9668849	1.0000000	NA	NA

pairwise <- pairwise.adonis(phylo_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.05690611	0.7893300	0.04999136	0.536	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.10866997	3.0547314	0.16919285	0.012	0.180
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.11760287	2.9988032	0.16661126	0.015	0.225
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09624422	1.7247115	0.10968160	0.175	1.000
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.18586788	4.1904227	0.21836011	0.005	0.075
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.23108846	4.0521838	0.20208192	0.004	0.060
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.26358465	4.3608960	0.21417997	0.004	0.060
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.25319427	3.2738422	0.17915456	0.038	0.570
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.39050120	5.9837393	0.27218933	0.001	0.015	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	0.14203376	5.4200212	0.25303529	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09666753	2.3682173	0.13635351	0.023	0.345
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.09252600	2.9824958	0.15711821	0.011	0.165
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.01842535	0.4144162	0.02688498	0.737	1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.05987967	1.7387847	0.09802164	0.115	1.000
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.03212966	0.6477782	0.04139746	0.694	1.000

6.3.4.1.5 Functional

func_post3 <- as.matrix(beta_q1f$S)
func_post3 <- as.dist(func_post3[rownames(func_post3) %in% samples_to_keep_post3,
               colnames(func_post3) %in% samples_to_keep_post3])
betadisper(func_post3, subset_meta_post3$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00125 0.0012477 0.0665    999  0.798
Residuals 50 0.93781 0.0187562                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.799
Hot_dry   0.79753        

adonis2(func_post3 ~ Population*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(func_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(func_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	0.03282552	0.02249363	0.3681798	0.363
Residual	48	1.42649950	0.97750637	NA	NA
Total	51	1.45932502	1.00000000	NA	NA

pairwise <- pairwise.adonis(func_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.022138495	0.64325180	0.041120082	0.406	1.00
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.017467416	0.54472723	0.035042572	0.440	1.00
Control.1_Acclimation vs Control.5_Post-FMT1	1	-0.007536319	-0.20662255	-0.013967233	0.906	1.00
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.128319067	3.48711706	0.199410632	0.068	1.00
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.051950542	1.26548063	0.077801613	0.293	1.00
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.001373886	0.07238159	0.004503476	0.640	1.00
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.002173211	0.09402475	0.005842215	0.632	1.00
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.059119657	2.62461793	0.148917721	0.145	1.00
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.010911935	0.39816986	0.024281360	0.484	1.00
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	-0.002303582	-0.11016709	-0.006933181	0.748	1.00
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.058908140	2.91987617	0.162940644	0.165	1.00
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.005904278	0.23427876	0.014431116	0.553	1.00
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.116146557	4.72479132	0.239535681	0.066	0.99
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.050009298	1.70482607	0.096291602	0.237	1.00
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.012358590	0.42381202	0.027477774	0.482	1.00

beta_richness_nmds_post3 <- richness_post3 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post3, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post3 <- neutral_post3 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post3, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post3 <- phylo_post3 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post3, by = join_by(sample == Tube_code))

beta_func_nmds_post3 <- func_post3 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post3, by = join_by(sample == Tube_code))

6.3.4.2 Comparison between Acclimation vs Post-FMT2

6.3.4.2.1 Number of samples used

[1] 53

6.3.4.2.2 Richness

richness_post4 <- as.matrix(beta_q0n$S)
richness_post4 <- as.dist(richness_post4[rownames(richness_post4) %in% samples_to_keep_post4,
               colnames(richness_post4) %in% samples_to_keep_post4])
betadisper(richness_post4, subset_meta_post4$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.07832 0.078322 11.371    999  0.001 ***
Residuals 51 0.35129 0.006888                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.001
Hot_dry  0.0014303        

adonis2(richness_post4 ~ Population*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	3.279021	0.1937289	3.924535	0.001
Residual	49	13.646799	0.8062711	NA	NA
Total	52	16.925820	1.0000000	NA	NA

#Arrange of metadata dataframe
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(richness_post4),]

pairwise <- pairwise.adonis(richness_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.3657243	1.123239	0.06966584	0.236	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.2605860	4.943717	0.24788342	0.001	0.015	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.8072604	2.940901	0.16392161	0.001	0.015	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.4817387	1.660775	0.09968176	0.023	0.345
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.1179704	3.885459	0.20573812	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.3606630	5.087152	0.24124415	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.9130048	3.195028	0.16645080	0.001	0.015	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.5959230	1.984036	0.11032208	0.002	0.030	.
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.2747787	4.275366	0.21086503	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6397330	2.913695	0.15405213	0.001	0.015	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.4575447	6.224524	0.28007456	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3276169	1.412318	0.08111028	0.037	0.555
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.6463814	2.560441	0.13795154	0.001	0.015	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.4796256	1.916520	0.10696943	0.003	0.045	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.1305044	4.268317	0.21059061	0.001	0.015	.

6.3.4.2.3 Neutral

neutral_post4 <- as.matrix(beta_q1n$S)
neutral_post4 <- as.dist(neutral_post4[rownames(neutral_post4) %in% samples_to_keep_post4,
               colnames(neutral_post4) %in% samples_to_keep_post4])
betadisper(neutral_post4, subset_meta_post4$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.07811 0.078108 9.6342    999  0.003 **
Residuals 51 0.41347 0.008107                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.004
Hot_dry  0.0031133        

adonis2(neutral_post4 ~ Population*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	3.823089	0.2357666	5.038847	0.001
Residual	49	12.392477	0.7642334	NA	NA
Total	52	16.215567	1.0000000	NA	NA

pairwise <- pairwise.adonis(neutral_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.2759858	0.9928976	0.06208366	0.488	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.4146464	6.5078325	0.30257965	0.001	0.015	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	1.1524353	4.9536068	0.24825621	0.001	0.015	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.4748999	1.9609749	0.11561688	0.005	0.075
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.3503168	5.4081420	0.26499923	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.6347704	6.8326887	0.29925029	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	1.3540292	5.3398081	0.25022756	0.001	0.015	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.6311089	2.4041625	0.13063146	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.6125755	5.9825981	0.27215155	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6202327	3.1519868	0.16457754	0.001	0.015	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.5701179	7.6327037	0.32297209	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3634438	1.7083388	0.09647087	0.042	0.630
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	1.0227481	4.6483346	0.22511910	0.001	0.015	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.5010202	2.2065321	0.12119453	0.002	0.030	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.3619424	5.7710313	0.26507845	0.001	0.015	.

6.3.4.2.4 Phylogenetic

phylo_post4 <- as.matrix(beta_q1p$S)
phylo_post4 <- as.dist(phylo_post4[rownames(phylo_post4) %in% samples_to_keep_post4,
               colnames(phylo_post4) %in% samples_to_keep_post4])
betadisper(phylo_post4, subset_meta_post4$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq     F N.Perm Pr(>F)  
Groups     1 0.03098 0.030984 2.885    999   0.08 .
Residuals 51 0.54772 0.010740                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.088
Hot_dry    0.0955        

adonis2(phylo_post4 ~ Population*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	0.6466392	0.2442162	5.277784	0.001
Residual	49	2.0011759	0.7557838	NA	NA
Total	52	2.6478151	1.0000000	NA	NA

pairwise <- pairwise.adonis(phylo_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.05690611	0.789330	0.04999136	0.510	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.10866997	3.054731	0.16919285	0.012	0.180
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.15591805	4.379209	0.22597458	0.007	0.105
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.07099742	1.879364	0.11134092	0.123	1.000
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.18682367	4.878754	0.24542555	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.23108846	4.052184	0.20208192	0.003	0.045	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.36496892	6.396667	0.28560797	0.001	0.015	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.22628210	3.829222	0.19311005	0.017	0.255
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.34830814	5.846334	0.26761166	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.10002871	4.383624	0.21505615	0.001	0.015	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	0.12577510	5.060129	0.24027055	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.06334378	2.499774	0.13512455	0.024	0.360
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.05927454	2.382025	0.12958449	0.024	0.360
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.06906280	2.722460	0.14541146	0.006	0.090
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.11081709	4.043656	0.20174244	0.002	0.030	.

6.3.4.2.5 Functional

func_post4 <- as.matrix(beta_q1f$S)
func_post4 <- as.dist(func_post4[rownames(func_post4) %in% samples_to_keep_post4,
               colnames(func_post4) %in% samples_to_keep_post4])
betadisper(func_post4, subset_meta_post4$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     1 0.00002 0.0000207 0.001    999  0.978
Residuals 51 1.04916 0.0205717                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.979
Hot_dry    0.9748        

adonis2(func_post4 ~ Population*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(func_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(func_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	0.02728193	0.0197108	0.3284164	0.569
Residual	49	1.35682891	0.9802892	NA	NA
Total	52	1.38411084	1.0000000	NA	NA

pairwise <- pairwise.adonis(func_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	2.213850e-02	0.64325180	0.0411200822	0.443	1
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.746742e-02	0.54472723	0.0350425723	0.427	1
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	3.674404e-02	1.08497265	0.0674525642	0.322	1
Control.1_Acclimation vs Control.6_Post-FMT2	1	2.983798e-02	0.93187248	0.0584910835	0.354	1
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.407466e-02	0.28346685	0.0185472873	0.517	1
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.373886e-03	0.07238159	0.0045034760	0.645	1
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	-4.929805e-04	-0.02385162	-0.0014929516	0.713	1
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	2.443617e-03	0.12903840	0.0080003779	0.576	1
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	6.968380e-03	0.19647180	0.0121305310	0.576	1
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	3.866607e-04	0.02093980	0.0013070267	0.673	1
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	-5.633463e-05	-0.00336651	-0.0002104511	0.681	1
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	2.011448e-03	0.06046867	0.0037650628	0.734	1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	-8.527330e-03	-0.46290555	-0.0297935723	0.849	1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	-1.648721e-03	-0.04717131	-0.0029569243	0.890	1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	4.367477e-03	0.13147026	0.0081499244	0.707	1

beta_richness_nmds_post4 <- richness_post4 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post4, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post4 <- neutral_post4 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post4, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post4 <- phylo_post4 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post4, by = join_by(sample == Tube_code))

beta_func_nmds_post4 <- func_post4 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post4, by = join_by(sample == Tube_code))

6.3.4.3 Comparison between Acclimation vs Post-FMT1 and Post-FMT2

6.3.4.3.1 Number of samples used

[1] 79

6.3.4.3.2 Richness

richness_post6 <- as.matrix(beta_q0n$S)
richness_post6 <- as.dist(richness_post6[rownames(richness_post6) %in% samples_to_keep_post6,
               colnames(richness_post6) %in% samples_to_keep_post6])
betadisper(richness_post6, subset_meta_post6$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.09792 0.097923 18.607    999  0.001 ***
Residuals 77 0.40523 0.005263                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.001
Hot_dry  4.714e-05        

adonis2(richness_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	8	6.395467	0.2451322	2.841433	0.001
Residual	70	19.694403	0.7548678	NA	NA
Total	78	26.089870	1.0000000	NA	NA

#Arrange of metadata dataframe
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(richness_post6),]

pairwise <- pairwise.adonis(richness_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.3657243	1.123239	0.06966584	0.218	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.2605860	4.943717	0.24788342	0.001	0.036	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.6804988	2.100611	0.12283837	0.002	0.072
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8481258	2.673290	0.16033371	0.001	0.036	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.1180661	3.809275	0.20252111	0.001	0.036	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.8072604	2.940901	0.16392161	0.001	0.036	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.4817387	1.660775	0.09968176	0.027	0.972
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.1179704	3.885459	0.20573812	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.3606630	5.087152	0.24124415	0.002	0.072
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.7216200	2.172734	0.11956009	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9551308	2.926054	0.16322910	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.2263345	4.039487	0.20157637	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.9130048	3.195028	0.16645080	0.001	0.036	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.5959230	1.984036	0.11032208	0.002	0.072
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.2747787	4.275366	0.21086503	0.001	0.036	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.4319792	5.384836	0.25180628	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8172413	3.194690	0.17558364	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.5796135	2.441615	0.13239702	0.002	0.072
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6397330	2.913695	0.15405213	0.001	0.036	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.4575447	6.224524	0.28007456	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3276169	1.412318	0.08111028	0.037	1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.5615418	1.729004	0.10335366	0.017	0.612
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.8438429	2.793772	0.14865413	0.002	0.072
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.7628135	2.683925	0.14364890	0.001	0.036	.
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.3432605	1.148733	0.06698647	0.246	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	1.1269580	3.799256	0.19188884	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.3734921	1.268929	0.07799710	0.115	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.3571397	1.297184	0.07959561	0.125	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.7769467	2.670898	0.15114670	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.6502360	2.253407	0.13060650	0.001	0.036	.
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.4132091	1.616138	0.09174188	0.009	0.324
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.0163992	3.760571	0.19030682	0.001	0.036	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.2732563	1.019281	0.05988979	0.429	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.6463814	2.560441	0.13795154	0.002	0.072
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.4796256	1.916520	0.10696943	0.004	0.144
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.1305044	4.268317	0.21059061	0.001	0.036	.

6.3.4.3.3 Neutral

neutral_post6 <- as.matrix(beta_q1n$S)
neutral_post6 <- as.dist(neutral_post6[rownames(neutral_post6) %in% samples_to_keep_post6,
               colnames(neutral_post6) %in% samples_to_keep_post6])
betadisper(neutral_post6, subset_meta_post6$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.09524 0.095237 15.396    999  0.001 ***
Residuals 77 0.47632 0.006186                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Cold_wet Hot_dry
Cold_wet              0.001
Hot_dry  0.00018818        

adonis2(neutral_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	8	7.736408	0.3121974	3.971673	0.001
Residual	70	17.044094	0.6878026	NA	NA
Total	78	24.780501	1.0000000	NA	NA

pairwise <- pairwise.adonis(neutral_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.2759858	0.9928976	0.06208366	0.482	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.4146464	6.5078325	0.30257965	0.001	0.036	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.8153894	3.0970603	0.17113610	0.001	0.036	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	1.1809241	4.4856470	0.24265567	0.002	0.072
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.4321524	5.7774260	0.27806264	0.001	0.036	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	1.1524353	4.9536068	0.24825621	0.001	0.036	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.4748999	1.9609749	0.11561688	0.008	0.288
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.3503168	5.4081420	0.26499923	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.6347704	6.8326887	0.29925029	0.001	0.036	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.9517634	3.3715700	0.17404733	0.002	0.072
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	1.3127773	4.6298256	0.23585668	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.6713369	6.2395460	0.28056085	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	1.3540292	5.3398081	0.25022756	0.001	0.036	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.6311089	2.4041625	0.13063146	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.6125755	5.9825981	0.27215155	0.001	0.036	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.5409781	6.8338056	0.29928456	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9133614	4.0964534	0.21451383	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.6954835	3.2951234	0.17077493	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6202327	3.1519868	0.16457754	0.001	0.036	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.5701179	7.6327037	0.32297209	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3634438	1.7083388	0.09647087	0.029	1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.6051778	2.2508491	0.13047758	0.016	0.576
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	1.0528902	4.1436369	0.20570451	0.001	0.036	.
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.8908158	3.7146920	0.18842252	0.003	0.108
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.3860927	1.5521758	0.08843210	0.081	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	1.3122237	5.1302726	0.24279254	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.4150076	1.6372683	0.09840968	0.047	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.3157079	1.3252026	0.08117526	0.169	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.0579520	4.2700097	0.22158835	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.7454015	2.9200493	0.16294873	0.002	0.072
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.4377161	1.9421261	0.10824392	0.010	0.360
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.3766597	5.8752789	0.26858075	0.001	0.036	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.3176516	1.3161367	0.07600637	0.147	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	1.0227481	4.6483346	0.22511910	0.001	0.036	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.5010202	2.2065321	0.12119453	0.003	0.108
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.3619424	5.7710313	0.26507845	0.001	0.036	.

6.3.4.3.4 Phylogenetic

phylo_post6 <- as.matrix(beta_q1p$S)
phylo_post6 <- as.dist(phylo_post6[rownames(phylo_post6) %in% samples_to_keep_post6,
               colnames(phylo_post6) %in% samples_to_keep_post6])
betadisper(phylo_post6, subset_meta_post6$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02554 0.025541 2.1633    999  0.165
Residuals 77 0.90909 0.011806                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.165
Hot_dry   0.14542        

adonis2(phylo_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	8	1.055866	0.2692287	3.223651	0.001
Residual	70	2.865952	0.7307713	NA	NA
Total	78	3.921818	1.0000000	NA	NA

pairwise <- pairwise.adonis(phylo_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.05690611	0.7893300	0.04999136	0.548	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.10866997	3.0547314	0.16919285	0.009	0.324
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.11760287	2.9988032	0.16661126	0.024	0.864
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09624422	1.7247115	0.10968160	0.165	1.000
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.18586788	4.1904227	0.21836011	0.005	0.180
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.15591805	4.3792085	0.22597458	0.007	0.252
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.07099742	1.8793639	0.11134092	0.122	1.000
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.18682367	4.8787543	0.24542555	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.23108846	4.0521838	0.20208192	0.002	0.072
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.26358465	4.3608960	0.21417997	0.009	0.324
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.25319427	3.2738422	0.17915456	0.048	1.000
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.39050120	5.9837393	0.27218933	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.36496892	6.3966666	0.28560797	0.001	0.036	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.22628210	3.8292220	0.19311005	0.020	0.720
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.34830814	5.8463335	0.26761166	0.001	0.036	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	0.14203376	5.4200212	0.25303529	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09666753	2.3682173	0.13635351	0.019	0.684
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.09252600	2.9824958	0.15711821	0.005	0.180
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.10002871	4.3836237	0.21505615	0.001	0.036	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	0.12577510	5.0601287	0.24027055	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.06334378	2.4997737	0.13512455	0.023	0.828
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.01842535	0.4144162	0.02688498	0.778	1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.05987967	1.7387847	0.09802164	0.120	1.000
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.07917244	3.0180046	0.15869197	0.009	0.324
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.04335491	1.5335604	0.08746429	0.189	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.10783045	3.7500438	0.18987521	0.002	0.072
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.03212966	0.6477782	0.04139746	0.682	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.06393539	1.5651817	0.09448624	0.138	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.05265949	1.2240203	0.07544494	0.287	1.000
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.09753501	2.2402429	0.12994265	0.015	0.540
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.07228545	2.3279593	0.12701683	0.037	1.000
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.11759094	3.5538444	0.18174658	0.001	0.036	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.06667255	1.9859527	0.11041687	0.096	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.05927454	2.3820253	0.12958449	0.024	0.864
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.06906280	2.7224602	0.14541146	0.004	0.144
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.11081709	4.0436561	0.20174244	0.002	0.072

6.3.4.3.5 Functional

func_post6 <- as.matrix(beta_q1f$S)
func_post6 <- as.dist(func_post6[rownames(func_post6) %in% samples_to_keep_post6,
               colnames(func_post6) %in% samples_to_keep_post6])
betadisper(func_post6, subset_meta_post6$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00012 0.0001211 0.0064    999  0.931
Residuals 77 1.45678 0.0189193                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.936
Hot_dry   0.93644        

adonis2(func_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(func_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(func_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	8	0.1638096	0.07637627	0.7235547	0.563
Residual	70	1.9809613	0.92362373	NA	NA
Total	78	2.1447709	1.00000000	NA	NA

pairwise <- pairwise.adonis(func_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	2.213850e-02	0.6432517992	4.112008e-02	0.444	1
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.746742e-02	0.5447272282	3.504257e-02	0.409	1
Control.1_Acclimation vs Control.5_Post-FMT1	1	-7.536319e-03	-0.2066225492	-1.396723e-02	0.886	1
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	1.283191e-01	3.4871170642	1.994106e-01	0.068	1
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	5.195054e-02	1.2654806328	7.780161e-02	0.291	1
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	3.674404e-02	1.0849726511	6.745256e-02	0.315	1
Control.1_Acclimation vs Control.6_Post-FMT2	1	2.983798e-02	0.9318724842	5.849108e-02	0.385	1
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.407466e-02	0.2834668511	1.854729e-02	0.497	1
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.373886e-03	0.0723815851	4.503476e-03	0.621	1
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	2.173211e-03	0.0940247450	5.842215e-03	0.600	1
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	5.911966e-02	2.6246179279	1.489177e-01	0.162	1
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.091194e-02	0.3981698590	2.428136e-02	0.501	1
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	-4.929805e-04	-0.0238516159	-1.492952e-03	0.714	1
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	2.443617e-03	0.1290384018	8.000378e-03	0.592	1
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	6.968380e-03	0.1964718029	1.213053e-02	0.547	1
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	-2.303582e-03	-0.1101670946	-6.933181e-03	0.747	1
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	5.890814e-02	2.9198761676	1.629406e-01	0.147	1
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	5.904278e-03	0.2342787606	1.443112e-02	0.548	1
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	3.866607e-04	0.0209397959	1.307027e-03	0.647	1
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	-5.633463e-05	-0.0033665098	-2.104511e-04	0.702	1
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	2.011448e-03	0.0604686736	3.765063e-03	0.741	1
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	1.161466e-01	4.7247913225	2.395357e-01	0.078	1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	5.000930e-02	1.7048260738	9.629160e-02	0.240	1
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	2.292148e-02	1.0143427269	5.961692e-02	0.324	1
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.305120e-02	0.6254809273	3.762182e-02	0.417	1
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	-8.269058e-03	-0.2211195871	-1.401364e-02	0.921	1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	1.235859e-02	0.4238120210	2.747777e-02	0.480	1
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	2.837514e-02	1.2912549771	7.926062e-02	0.288	1
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	3.882231e-02	1.9287804257	1.139350e-01	0.227	1
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	5.671762e-02	1.5020399981	9.102147e-02	0.256	1
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	-3.766327e-03	-0.1400672850	-8.831518e-03	0.725	1
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	2.020599e-03	0.0803166527	4.994718e-03	0.618	1
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	5.871214e-06	0.0001408358	8.802161e-06	0.857	1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	-8.527330e-03	-0.4629055475	-2.979357e-02	0.875	1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	-1.648721e-03	-0.0471713074	-2.956924e-03	0.900	1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	4.367477e-03	0.1314702630	8.149924e-03	0.710	1

beta_richness_nmds_post6 <- richness_post6 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post6, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post6 <- neutral_post6 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post6, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post6 <- phylo_post6 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post6, by = join_by(sample == Tube_code))

beta_func_nmds_post6 <- func_post6 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post6, by = join_by(sample == Tube_code))

6.3.5 5. Are there differences between the control and the treatment group?

6.3.5.1 After 1 week –> Post-FMT1

6.3.5.1.1 Number of samples used

[1] 26

6.3.5.1.2 Richness

richness_post1 <- as.matrix(beta_q0n$S)
richness_post1 <- as.dist(richness_post1[rownames(richness_post1) %in% samples_to_keep_post1,
               colnames(richness_post1) %in% samples_to_keep_post1])
betadisper(richness_post1, subset_meta_post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.017675 0.0088373 2.3825    999  0.106
Residuals 23 0.085312 0.0037092                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0050000     0.652
Hot_control 0.0068795                 0.196
Treatment   0.6248469   0.2084296          

adonis2(richness_post1 ~ type,
        data = subset_meta_post1 %>% arrange(match(Tube_code,labels(richness_post1))),
        permutations = 999,
        strata = subset_meta_post1 %>% arrange(match(Tube_code,labels(richness_post1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.195567	0.1448246	1.947534	1
Residual	23	7.059710	0.8551754	NA	NA
Total	25	8.255277	1.0000000	NA	NA

#Arrange of metadata dataframe
subset_meta_post1_arrange <- column_to_rownames(subset_meta_post1, "Tube_code")
subset_meta_post1_arrange<-subset_meta_post1_arrange[labels(richness_post1),]

pairwise <- pairwise.adonis(richness_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control vs Treatment	1	0.5615418	1.729004	0.1033537	0.013	0.039	.
Control vs Hot_control	1	0.8438429	2.793772	0.1486541	0.001	0.003	*
Treatment vs Hot_control	1	0.3734921	1.268929	0.0779971	0.116	0.348

6.3.5.1.3 Neutral

neutral_post1 <- as.matrix(beta_q1n$S)
neutral_post1 <- as.dist(neutral_post1[rownames(neutral_post1) %in% samples_to_keep_post1,
               colnames(neutral_post1) %in% samples_to_keep_post1])
betadisper(neutral_post1, subset_meta_post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.011001 0.0055005 0.6303    999  0.563
Residuals 23 0.200714 0.0087267                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.22600     0.945
Hot_control 0.21166                 0.452
Treatment   0.95468     0.43604          

adonis2(neutral_post1 ~ type,
        data = subset_meta_post1 %>% arrange(match(Tube_code,labels(neutral_post1))),
        permutations = 999,
        strata = subset_meta_post1 %>% arrange(match(Tube_code,labels(neutral_post1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.395968	0.1900228	2.697931	1
Residual	23	5.950350	0.8099772	NA	NA
Total	25	7.346318	1.0000000	NA	NA

pairwise <- pairwise.adonis(neutral_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control vs Treatment	1	0.6051778	2.250849	0.13047758	0.012	0.036	.
Control vs Hot_control	1	1.0528902	4.143637	0.20570451	0.001	0.003	*
Treatment vs Hot_control	1	0.4150076	1.637268	0.09840968	0.053	0.159

6.3.5.1.4 Phylogenetic

phylo_post1 <- as.matrix(beta_q1p$S)
phylo_post1 <- as.dist(phylo_post1[rownames(phylo_post1) %in% samples_to_keep_post1,
               colnames(phylo_post1) %in% samples_to_keep_post1])
betadisper(phylo_post1, subset_meta_post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00440 0.0021994 0.1369    999  0.898
Residuals 23 0.36941 0.0160614                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.91500     0.690
Hot_control 0.91505                 0.788
Treatment   0.63312     0.73046          

adonis2(phylo_post1 ~ type,
        data = subset_meta_post1 %>% arrange(match(Tube_code,labels(phylo_post1))),
        permutations = 999,
        strata = subset_meta_post1 %>% arrange(match(Tube_code,labels(phylo_post1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	0.07451036	0.07059466	0.8735033	1
Residual	23	0.98095695	0.92940534	NA	NA
Total	25	1.05546731	1.00000000	NA	NA

pairwise <- pairwise.adonis(phylo_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control vs Treatment	1	0.01842535	0.4144162	0.02688498	0.792	1.00
Control vs Hot_control	1	0.05987967	1.7387847	0.09802164	0.130	0.39
Treatment vs Hot_control	1	0.03212966	0.6477782	0.04139746	0.701	1.00

6.3.5.1.5 Functional

func_post1 <- as.matrix(beta_q1f$S)
func_post1 <- as.dist(func_post1[rownames(func_post1) %in% samples_to_keep_post1,
               colnames(func_post1) %in% samples_to_keep_post1])
betadisper(func_post1, subset_meta_post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     2 0.00400 0.0020014 0.145    999  0.862
Residuals 23 0.31753 0.0138057                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.59500     0.732
Hot_control 0.59817                 0.821
Treatment   0.75141     0.83718          

adonis2(func_post1 ~ type,
        data = subset_meta_post1 %>% arrange(match(Tube_code,labels(func_post1))),
        permutations = 999,
        strata = subset_meta_post1 %>% arrange(match(Tube_code,labels(func_post1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	0.1186444	0.1568566	2.139435	1
Residual	23	0.6377435	0.8431434	NA	NA
Total	25	0.7563879	1.0000000	NA	NA

pairwise <- pairwise.adonis(func_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control vs Treatment	1	0.11614656	4.724791	0.23953568	0.072	0.216
Control vs Hot_control	1	0.05000930	1.704826	0.09629160	0.228	0.684
Treatment vs Hot_control	1	0.01235859	0.423812	0.02747777	0.507	1.000

beta_richness_nmds_post1 <- richness_post1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post1, by = join_by(sample == Tube_code))

beta_neutral_nmds_post1 <- neutral_post1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post1, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post1 <- phylo_post1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post1, by = join_by(sample == Tube_code))

beta_functional_nmds_post1 <- func_post1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post1, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.5.2 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")

6.3.5.2.1 Number of samples used

[1] 27

6.3.5.2.2 Richness

richness_post2 <- as.matrix(beta_q0n$S)
richness_post2 <- as.dist(richness_post2[rownames(richness_post2) %in% samples_to_keep_post2,
               colnames(richness_post2) %in% samples_to_keep_post2])
betadisper(richness_post2, subset_meta_post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.002011 0.0010056 0.1982    999  0.845
Residuals 24 0.121775 0.0050740                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.70000     0.827
Hot_control 0.67789                 0.633
Treatment   0.79246     0.59820          

adonis2(richness_post2 ~ type,
        data = subset_meta_post2 %>% arrange(match(Tube_code,labels(richness_post2))),
        permutations = 999,
        strata = subset_meta_post2 %>% arrange(match(Tube_code,labels(richness_post2))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.504341	0.1967776	2.939822	1
Residual	24	6.140538	0.8032224	NA	NA
Total	26	7.644879	1.0000000	NA	NA

#Arrange of metadata dataframe
subset_meta_post2_arrange <- column_to_rownames(subset_meta_post2, "Tube_code")
subset_meta_post2_arrange<-subset_meta_post2_arrange[labels(richness_post2),]

pairwise <- pairwise.adonis(richness_post2, subset_meta_post2_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	0.6463814	2.560441	0.1379515	0.001	0.003	*
Treatment vs Hot_control	1	0.4796256	1.916520	0.1069694	0.001	0.003	*
Control vs Hot_control	1	1.1305044	4.268317	0.2105906	0.001	0.003	*

6.3.5.2.3 Neutral

neutral_post2 <- as.matrix(beta_q1n$S)
neutral_post2 <- as.dist(neutral_post2[rownames(neutral_post2) %in% samples_to_keep_post2,
               colnames(neutral_post2) %in% samples_to_keep_post2])
betadisper(neutral_post2, subset_meta_post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.008262 0.0041311 0.8024    999  0.458
Residuals 24 0.123559 0.0051483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.44000     0.649
Hot_control 0.44675                 0.249
Treatment   0.65989     0.25095          

adonis2(neutral_post2 ~ type,
        data = subset_meta_post2 %>% arrange(match(Tube_code,labels(neutral_post2))),
        permutations = 999,
        strata = subset_meta_post2 %>% arrange(match(Tube_code,labels(neutral_post2))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.923807	0.2603795	4.224537	1
Residual	24	5.464666	0.7396205	NA	NA
Total	26	7.388473	1.0000000	NA	NA

pairwise <- pairwise.adonis(neutral_post2, subset_meta_post2_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	1.0227481	4.648335	0.2251191	0.001	0.003	*
Treatment vs Hot_control	1	0.5010202	2.206532	0.1211945	0.003	0.009	*
Control vs Hot_control	1	1.3619424	5.771031	0.2650785	0.001	0.003	*

6.3.5.2.4 Phylogenetic

phylo_post2 <- as.matrix(beta_q1p$S)
phylo_post2 <- as.dist(phylo_post2[rownames(phylo_post2) %in% samples_to_keep_post2,
               colnames(phylo_post2) %in% samples_to_keep_post2])
betadisper(phylo_post2, subset_meta_post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.000407 0.0002034 0.0487    999  0.947
Residuals 24 0.100305 0.0041794                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.92800     0.846
Hot_control 0.93765                 0.737
Treatment   0.83933     0.76015          

adonis2(phylo_post2 ~ type,
        data = subset_meta_post2 %>% arrange(match(Tube_code,labels(phylo_post2))),
        permutations = 999,
        strata = subset_meta_post2 %>% arrange(match(Tube_code,labels(phylo_post2))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	0.1594363	0.2042241	3.079623	1
Residual	24	0.6212564	0.7957759	NA	NA
Total	26	0.7806927	1.0000000	NA	NA

pairwise <- pairwise.adonis(phylo_post2, subset_meta_post2_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	0.05927454	2.382025	0.1295845	0.027	0.081
Treatment vs Hot_control	1	0.06906280	2.722460	0.1454115	0.004	0.012	.
Control vs Hot_control	1	0.11081709	4.043656	0.2017424	0.002	0.006	*

6.3.5.2.5 Functional

func_post2 <- as.matrix(beta_q1f$S)
func_post2 <- as.dist(func_post2[rownames(func_post2) %in% samples_to_keep_post2,
               colnames(func_post2) %in% samples_to_keep_post2])
betadisper(func_post2, subset_meta_post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01126 0.0056302 0.2861    999  0.813
Residuals 24 0.47233 0.0196806                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.55500     0.681
Hot_control 0.48255                 0.816
Treatment   0.60116     0.75643          

adonis2(func_post2 ~ type,
        data = subset_meta_post2 %>% arrange(match(Tube_code,labels(func_post2))),
        permutations = 999,
        strata = subset_meta_post2 %>% arrange(match(Tube_code,labels(func_post2))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	-0.003872383	-0.005621319	-0.06707875	1
Residual	24	0.692746821	1.005621319	NA	NA
Total	26	0.688874438	1.000000000	NA	NA

pairwise <- pairwise.adonis(func_post2, subset_meta_post2_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	-0.008527330	-0.46290555	-0.029793572	0.836	1
Treatment vs Hot_control	1	-0.001648721	-0.04717131	-0.002956924	0.906	1
Control vs Hot_control	1	0.004367477	0.13147026	0.008149924	0.699	1

beta_richness_nmds_post2 <- richness_post2 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post2, by = join_by(sample == Tube_code))

beta_neutral_nmds_post2 <- neutral_post2 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post2, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post2 <- phylo_post2 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post2, by = join_by(sample == Tube_code))

beta_functional_nmds_post2 <- func_post2 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post2, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.5.3 Post1 vs Post2

post5 <- meta %>%
  filter(time_point == "6_Post-FMT2" | time_point == "5_Post-FMT1")

post5.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post5))]
identical(sort(colnames(post5.counts)),sort(as.character(rownames(post5))))

post5_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2"| time_point == "5_Post-FMT1")

6.3.5.3.1 Number of samples used

[1] 53

6.3.5.3.2 Richness

richness_post5 <- as.matrix(beta_q0n$S)
richness_post5 <- as.dist(richness_post5[rownames(richness_post5) %in% samples_to_keep_post5,
               colnames(richness_post5) %in% samples_to_keep_post5])
betadisper(richness_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01841 0.0092048 1.7364    999  0.178
Residuals 50 0.26505 0.0053010                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
             Control Hot_control Treatment
Control                 0.033000     0.714
Hot_control 0.039117                 0.229
Treatment   0.716358    0.218648          

adonis2(richness_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	3.261916	0.1981463	2.322836	0.001
Residual	47	13.200248	0.8018537	NA	NA
Total	52	16.462164	1.0000000	NA	NA

#Arrange of metadata dataframe
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(richness_post5),]

pairwise <- pairwise.adonis(richness_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.5615418	1.729004	0.10335366	0.020	0.300
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.8438429	2.793772	0.14865413	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.7628135	2.683925	0.14364890	0.001	0.015	.
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.3432605	1.148733	0.06698647	0.274	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	1.1269580	3.799256	0.19188884	0.002	0.030	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.3734921	1.268929	0.07799710	0.117	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.3571397	1.297184	0.07959561	0.128	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.7769467	2.670898	0.15114670	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.6502360	2.253407	0.13060650	0.001	0.015	.
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.4132091	1.616138	0.09174188	0.008	0.120
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.0163992	3.760571	0.19030682	0.001	0.015	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.2732563	1.019281	0.05988979	0.413	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.6463814	2.560441	0.13795154	0.001	0.015	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.4796256	1.916520	0.10696943	0.001	0.015	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.1305044	4.268317	0.21059061	0.001	0.015	.

6.3.5.3.3 Neutral

neutral_post5 <- as.matrix(beta_q1n$S)
neutral_post5 <- as.dist(neutral_post5[rownames(neutral_post5) %in% samples_to_keep_post5,
               colnames(neutral_post5) %in% samples_to_keep_post5])
betadisper(neutral_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     2 0.01992 0.0099587 1.565    999  0.202
Residuals 50 0.31818 0.0063636                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.10400     0.881
Hot_control 0.10701                 0.180
Treatment   0.87156     0.17449          

adonis2(neutral_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	3.947979	0.2569798	3.251069	0.001
Residual	47	11.415016	0.7430202	NA	NA
Total	52	15.362995	1.0000000	NA	NA

pairwise <- pairwise.adonis(neutral_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.6051778	2.250849	0.13047758	0.015	0.225
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	1.0528902	4.143637	0.20570451	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.8908158	3.714692	0.18842252	0.002	0.030	.
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.3860927	1.552176	0.08843210	0.077	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	1.3122237	5.130273	0.24279254	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.4150076	1.637268	0.09840968	0.054	0.810
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.3157079	1.325203	0.08117526	0.158	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.0579520	4.270010	0.22158835	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.7454015	2.920049	0.16294873	0.001	0.015	.
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.4377161	1.942126	0.10824392	0.008	0.120
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.3766597	5.875279	0.26858075	0.001	0.015	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.3176516	1.316137	0.07600637	0.182	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	1.0227481	4.648335	0.22511910	0.001	0.015	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.5010202	2.206532	0.12119453	0.002	0.030	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.3619424	5.771031	0.26507845	0.001	0.015	.

6.3.5.3.4 Phylogenetic

phylo_post5 <- as.matrix(beta_q1p$S)
phylo_post5 <- as.dist(phylo_post5[rownames(phylo_post5) %in% samples_to_keep_post5,
               colnames(phylo_post5) %in% samples_to_keep_post5])
betadisper(phylo_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00051 0.0002543 0.0265    999   0.97
Residuals 50 0.47996 0.0095993                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.87900     0.818
Hot_control 0.88926                 0.928
Treatment   0.82391     0.91902          

adonis2(phylo_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	0.3518222	0.180049	2.0641	0.001
Residual	47	1.6022134	0.819951	NA	NA
Total	52	1.9540356	1.000000	NA	NA

pairwise <- pairwise.adonis(phylo_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.01842535	0.4144162	0.02688498	0.806	1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.05987967	1.7387847	0.09802164	0.117	1.000
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.07917244	3.0180046	0.15869197	0.007	0.105
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.04335491	1.5335604	0.08746429	0.162	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.10783045	3.7500438	0.18987521	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.03212966	0.6477782	0.04139746	0.669	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.06393539	1.5651817	0.09448624	0.126	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.05265949	1.2240203	0.07544494	0.303	1.000
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.09753501	2.2402429	0.12994265	0.010	0.150
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.07228545	2.3279593	0.12701683	0.040	0.600
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.11759094	3.5538444	0.18174658	0.001	0.015	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.06667255	1.9859527	0.11041687	0.090	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.05927454	2.3820253	0.12958449	0.022	0.330
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.06906280	2.7224602	0.14541146	0.004	0.060
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.11081709	4.0436561	0.20174244	0.002	0.030	.

6.3.5.3.5 Functional

func_post5 <- as.matrix(beta_q1f$S)
func_post5 <- as.dist(func_post5[rownames(func_post5) %in% samples_to_keep_post5,
               colnames(func_post5) %in% samples_to_keep_post5])
betadisper(func_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00855 0.0042775 0.2522    999  0.787
Residuals 50 0.84815 0.0169631                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.52500     0.611
Hot_control 0.50303                 0.844
Treatment   0.59083     0.82686          

adonis2(func_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(func_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(func_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	0.1047944	0.07301299	0.7403794	0.471
Residual	47	1.3304903	0.92698701	NA	NA
Total	52	1.4352848	1.00000000	NA	NA

pairwise <- pairwise.adonis(func_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	1.161466e-01	4.7247913225	2.395357e-01	0.075	1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	5.000930e-02	1.7048260738	9.629160e-02	0.230	1
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	2.292148e-02	1.0143427269	5.961692e-02	0.311	1
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.305120e-02	0.6254809273	3.762182e-02	0.404	1
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	-8.269058e-03	-0.2211195871	-1.401364e-02	0.921	1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	1.235859e-02	0.4238120210	2.747777e-02	0.508	1
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	2.837514e-02	1.2912549771	7.926062e-02	0.301	1
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	3.882231e-02	1.9287804257	1.139350e-01	0.199	1
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	5.671762e-02	1.5020399981	9.102147e-02	0.242	1
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	-3.766327e-03	-0.1400672850	-8.831518e-03	0.731	1
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	2.020599e-03	0.0803166527	4.994718e-03	0.607	1
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	5.871214e-06	0.0001408358	8.802161e-06	0.853	1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	-8.527330e-03	-0.4629055475	-2.979357e-02	0.859	1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	-1.648721e-03	-0.0471713074	-2.956924e-03	0.910	1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	4.367477e-03	0.1314702630	8.149924e-03	0.696	1

beta_richness_nmds_post5 <- richness_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))

beta_neutral_nmds_post5 <- neutral_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post5 <- phylo_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))

beta_functional_nmds_post5 <- func_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")